OHI British Columbia | OHI Science | Citation policy
Currently, the Lasting Special Places goal model is identical to the OHI Global model: a region’s status is based upon percent of protected area within 1 km inland buffer and percent of protected area within 3 nautical mile offshore buffer, compared to a reference point of 30% protected area.
\[X_{LSP} = \frac{\frac{pA_{CMPA}}{pA_{refCMPA}} + \frac{pA_{CP}}{pA_{refCP}}}{2}\]
pA = percent of area within the inland or offshore buffer; CMPA = coastal marine protected area (3nm offshore); CP = coastline protected (1km inland); and refCMPA = refCP = 30% reference point for both measures.
Future changes may incorporate other data sets and MaPP planning zones.
WDPA data base
If the WDPA-MPA shapefile has already been rasterized for the global assessment, then we may be able to use that directly. It is Mollweide projection which is not ideal for BC, as we would prefer to work in BC Albers. Some potential solutions:
We’re going with the fourth option.
NOTE: If BC WDPA file does not yet exist, get_wdpa_poly() creates it from the original WDPA-MPA file. This takes a long time, due to reading in the full WDPA-MPA geodatabase into a SpatialPolygonsDataFrame.
poly_wdpa_bc <- get_wdpa_poly(p4s_bcalb, reload = FALSE) ### defaults to BC Albersrast_eez <- raster(file.path(dir_spatial, 'ohibc_rgn_raster_500m.tif'))
wdpa_bc_shp_file <- file.path(dir_goal_anx, 'int', 'wdpa_bc_bcalb.shp')
wdpa_bc_rast_file <- file.path(dir_rast, 'rast_wdpa_bc.tif')
lsp_rasterize(wdpa_bc_shp_file,
wdpa_bc_rast_file,
rast_eez, 'STATUS_YR')## class : RasterBrick
## dimensions : 3126, 3424, 10703424, 1 (nrow, ncol, ncell, nlayers)
## resolution : 500, 500 (x, y)
## extent : 159000, 1871000, 173000, 1736000 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## data source : /home/ohara/github/ohibc/prep/lsp/v2016/raster/rast_wdpa_bc_tmp.tif
## names : rast_wdpa_bc_tmp
## NULL
rast_wdpa_bc <- raster::raster(wdpa_bc_rast_file)Buffer shapefiles are located in github/ohibc/prep/spatial. LSP uses 1 km inland and 3nm offshore buffers, while resilience requires analysis over the entire EEZ.
Analysis will be done using raster::crosstab() comparing the WDPA raster to various region rasters. Using a 500 m raster is the coarsest that should be used on a 1 km feature; a base raster is available at ~/github/ohibc/prep/spatial/ohibc_rgn_raster_500m.tif.
library(tmap)
rast_map <- tm_shape(rast_3nm, is.master = TRUE) +
tm_raster(alpha = 1, palette = 'Blues') +
tm_shape(rast_1km) +
tm_raster(alpha = 1, palette = 'Greens') +
tm_shape(rast_wdpa_bc) +
tm_raster(alpha = .5, palette = 'Reds')
print(rast_map)Once the WDPA raster is cross-tabulated against the OHI region rasters (both 3 nm offshore and 1 km inland) we have the number of protected cells, identified by year of protection, within each region. NA values are unprotected cells.
## year rgn_id n_cells
## Min. :1886 Min. :1.000 Min. : 0
## 1st Qu.:1947 1st Qu.:2.000 1st Qu.: 0
## Median :1970 Median :4.000 Median : 0
## Mean :1968 Mean :4.143 Mean : 15204
## 3rd Qu.:1992 3rd Qu.:6.000 3rd Qu.: 0
## Max. :2013 Max. :8.000 Max. :9942436
## NA's :8 NA's :88
## year rgn_id n_cells
## Min. :1886 Min. :1.000 Min. : 0
## 1st Qu.:1947 1st Qu.:2.000 1st Qu.: 0
## Median :1970 Median :4.000 Median : 0
## Mean :1968 Mean :4.143 Mean : 15204
## 3rd Qu.:1992 3rd Qu.:6.000 3rd Qu.: 3
## Max. :2013 Max. :8.000 Max. :10012355
## NA's :8 NA's :88
Grouping by rgn_id, the total number of cells per region is determined by summing cell counts across ALL years, including cells with year == NA (unprotected cells). We can then determine the protected area for each year by looking at the cumulative sum of cells up to any given year.
Since the cells are 500 m on a side, we can easily calculate area by multiplying cell count * 0.25 km2 per cell.
Finally we can calculate the status of a region for any given year by finding the ratio of protected:total and normalizing by the goal’s target of 30% protected area.
The status is based on a simple arithmetic average of the inland and offshore status values.
We can save outputs for the following layers:
a_prot_inland_file <- file.path(dir_goal, ‘output’, ‘lsp_protected_inland1km.csv’) a_prot_offshore_file <- file.path(dir_goal, ‘output’, ‘lsp_protected_offshore3nm.csv’) a_tot_inland_file <- file.path(dir_goal, ‘output’, ‘lsp_a_total_inland1km.csv’) a_tot_offshore_file <- file.path(dir_goal, ‘output’, ‘lsp_a_total_offshore3nm.csv’)
~/github/ohibc/prep/lsp/v2016/output/lsp_protected_inland1km.csv: inland protected area (km2) for each region (since 1980)~/github/ohibc/prep/lsp/v2016/output/lsp_protected_offshore3nm.csv: offshore protected area (km2) for each region (since 1980)~/github/ohibc/prep/lsp/v2016/output/lsp_a_total_inland1km.csv: inland 1 km total area (km2) for each region~/github/ohibc/prep/lsp/v2016/output/lsp_a_total_offshore3nm.csv: offshore 3 nm total area (km2) for each regionFrom these layers, we can also estimate the status and trend. “Official” values will be determined in the toolbox? Trend is based on linear model going back ten years from each status year to smooth trend values, since addition of new MPAs is rather sporadic.
Year-by-year status and trend estimates will be saved:
~/github/ohibc/prep/lsp/v2016/output/lsp_status.csv: estimate of status by region since 1980~/github/ohibc/prep/lsp/v2016/output/lsp_trend.csv: estimate of trend by region since 1990Examining OHIBC Lasting Special Places scores for 1995, 2005, and 2015:
To examine results, we plot the estimated status and trend over time.